
       m_a   = (*fields_ptr).u[M_A  ]GRD_IND;
       m_a_t = (*fields_ptr).u[M_A_T]GRD_IND;
       m_a_x = (*fields_ptr).u[M_A_X]GRD_IND;
       m_a_y = (*fields_ptr).u[M_A_Y]GRD_IND;

       m_b   = (*fields_ptr).u[M_B  ]GRD_IND;
       m_b_t = (*fields_ptr).u[M_B_T]GRD_IND;
       m_b_x = (*fields_ptr).u[M_B_X]GRD_IND;
       m_b_y = (*fields_ptr).u[M_B_Y]GRD_IND;

       m_c   = (*fields_ptr).u[M_C  ]GRD_IND;
       m_c_t = (*fields_ptr).u[M_C_T]GRD_IND;
       m_c_x = (*fields_ptr).u[M_C_X]GRD_IND;
       m_c_y = (*fields_ptr).u[M_C_Y]GRD_IND;

		// Metric components 
		
		g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
		g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;
       	g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;

       	Sdet = 1./sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);
       	
       	// first the wave equation
       	
       m_a_t_t =   (dfields.du[D_X_M_A_X]GRD_IND + dfields.du[D_Y_M_A_Y]GRD_IND)/Sdet;

       m_a_x_t = dfields.du[D_X_M_A_T]GRD_IND;
       m_a_y_t = dfields.du[D_Y_M_A_T]GRD_IND;

       m_b_t_t =   (dfields.du[D_X_M_B_X]GRD_IND + dfields.du[D_Y_M_B_Y]GRD_IND)/Sdet;

       m_b_x_t = dfields.du[D_X_M_B_T]GRD_IND;
       m_b_y_t = dfields.du[D_Y_M_B_T]GRD_IND;

       m_c_t_t =   (dfields.du[D_X_M_C_X]GRD_IND + dfields.du[D_Y_M_C_Y]GRD_IND)/Sdet;

       m_c_x_t = dfields.du[D_X_M_C_T]GRD_IND;
       m_c_y_t = dfields.du[D_Y_M_C_T]GRD_IND;



			norm_m_2 = (m_a*m_a + m_b*m_b + m_c*m_c);
			norm_m = sqrt(norm_m_2);


#ifdef STANDARD

#else

#ifdef FORCED_NORMALIZATION

#else
// Normalizing m_a

         m_a = m_a/norm_m;
         m_b = m_b/norm_m;
         m_c = m_c/norm_m;
 
 	(*fields_ptr).u[M_A]GRD_IND = m_a;
 	(*fields_ptr).u[M_B]GRD_IND = m_b;
 	(*fields_ptr).u[M_C]GRD_IND = m_c;

#endif

//// Normalizing m_a_t

         //m_t_dot_m = m_a*m_a_t + m_b*m_b_t + m_c*m_c_t;
 
         //m_a_t = m_a_t - m_a*m_t_dot_m;
         //m_b_t = m_b_t - m_b*m_t_dot_m;
         //m_c_t = m_c_t - m_c*m_t_dot_m;


//// Normalizing m_a_x

         //m_t_dot_m = m_a*m_a_x + m_b*m_b_x + m_c*m_c_x;
 
         //m_a_x = m_a_x - m_a*m_t_dot_m;
         //m_b_x = m_b_x - m_b*m_t_dot_m;
         //m_c_x = m_c_x - m_c*m_t_dot_m;

//// Normalizing m_a_y

         //m_t_dot_m = m_a*m_a_y + m_b*m_b_y + m_c*m_c_y;
 
         //m_a_y = m_a_y - m_a*m_t_dot_m;
         //m_b_y = m_b_y - m_b*m_t_dot_m;
         //m_c_y = m_c_y - m_c*m_t_dot_m;
 
#endif // STANDARD

#ifdef STANDARD

      d2_aa = (1.-m_a*m_a)*(g_up_xx*m_a_x*m_a_x +
                          2.*g_up_xy*m_a_x*m_a_y +
                             g_up_yy*m_a_y*m_a_y - m_a_t*m_a_t);

       d2_bb = (1.-m_b*m_b)*(g_up_xx*m_b_x*m_b_x +
                          2.*g_up_xy*m_b_x*m_b_y +
                             g_up_yy*m_b_y*m_b_y - m_b_t*m_b_t);

       d2_cc = (1.-m_c*m_c)*(g_up_xx*m_c_x*m_c_x +
                          2.*g_up_xy*m_c_x*m_c_y +
                             g_up_yy*m_c_y*m_c_y - m_c_t*m_c_t);

       d2_ab = -2.*m_a*m_b*(g_up_xx*m_a_x*m_b_x +
                         2.*g_up_xy*m_a_x*m_b_y +
                            g_up_yy*m_a_y*m_b_y - m_a_t*m_b_t);


       d2_ac = -2.*m_a*m_c*(g_up_xx*m_a_x*m_c_x +
                         2.*g_up_xy*m_a_x*m_c_y +
                            g_up_yy*m_a_y*m_c_y - m_a_t*m_c_t);


       d2_cb = -2.*m_c*m_b*(g_up_xx*m_c_x*m_b_x +
                         2.*g_up_xy*m_c_x*m_b_y +
                            g_up_yy*m_c_y*m_b_y - m_c_t*m_b_t);



#endif
    

#ifdef PROJECTION

       d2_aa = (g_up_xx*m_a_x*m_a_x +
                          2.*g_up_xy*m_a_x*m_a_y +
                             g_up_yy*m_a_y*m_a_y - m_a_t*m_a_t);

       d2_bb = (g_up_xx*m_b_x*m_b_x +
                          2.*g_up_xy*m_b_x*m_b_y +
                             g_up_yy*m_b_y*m_b_y - m_b_t*m_b_t);

       d2_cc = (g_up_xx*m_c_x*m_c_x +
                          2.*g_up_xy*m_c_x*m_c_y +
                             g_up_yy*m_c_y*m_c_y - m_c_t*m_c_t);

       d2_ab = 0.0;


       d2_ac = 0.0;


       d2_cb = 0.0;


#endif





//printf("ee=%f, mm=%f\n",ee,mm);

#ifdef PROJECTION

       d2 = -ee*(d2_aa + d2_bb + d2_cc + d2_ab + d2_ac + d2_cb);

       dip = 0.0; 
       dip_t = 0.0;
#endif

#ifdef STANDARD

       d2 = -ee*(d2_aa + d2_bb + d2_cc + d2_ab + d2_ac + d2_cb)/norm_m_2;
       
       dip =  -mm*(1. -norm_m_2);
       dip_t = mm*(1.-norm_m_2)*(1.-norm_m_2); 

#endif

#ifdef FULL_PROJECTION



// Normalizing m_a_t_t

         m_t_dot_m = (m_a*m_a_t_t + m_b*m_b_t_t + m_c*m_c_t_t) 
                   + (m_a_t*m_a_t + m_b_t*m_b_t + m_c_t*m_c_t) ;
 
         m_a_t_t = m_a_t_t - m_a*m_t_dot_m;
         m_b_t_t = m_b_t_t - m_b*m_t_dot_m;
         m_c_t_t = m_c_t_t - m_c*m_t_dot_m;
//printf("inside FP");


#else

       d2 = d2 + dip;

       m_a_t_t =  m_a_t_t - m_a*d2 - m_a_t*dip_t;

       m_b_t_t =  m_b_t_t - m_b*d2 - m_b_t*dip_t;

       m_c_t_t =  m_c_t_t - m_c*d2 - m_c_t*dip_t;

#endif


#ifdef DISSIPATION

  m_a_t_t = m_a_t_t + sigma*dfields.du[DISS_M_A_T]GRD_IND;
  m_a_x_t = m_a_x_t + sigma*dfields.du[DISS_M_A_X]GRD_IND;
  m_a_y_t = m_a_y_t + sigma*dfields.du[DISS_M_A_Y]GRD_IND;

  m_b_t_t = m_b_t_t + sigma*dfields.du[DISS_M_B_T]GRD_IND;
  m_b_x_t = m_b_x_t + sigma*dfields.du[DISS_M_B_X]GRD_IND;
  m_b_y_t = m_b_y_t + sigma*dfields.du[DISS_M_B_Y]GRD_IND;

  m_c_t_t = m_c_t_t + sigma*dfields.du[DISS_M_C_T]GRD_IND;
  m_c_x_t = m_c_x_t + sigma*dfields.du[DISS_M_C_X]GRD_IND;
  m_c_y_t = m_c_y_t + sigma*dfields.du[DISS_M_C_Y]GRD_IND;



#endif

/* printf("I am inside eq_eq"); */
